#!/usr/bin/perl
# 
# Exercise 10.7
# 
# Extend the program in Exercise 10.6 to count all the sequences in a GenBank 
# library.
#
# Answer to Exercise 10.7
#
# This is a straightforward extension of the use of the approach of Exercise 10.6
# with the examples of Chapter 10.  The necessary subroutines are in the module
# BeginPerlBioinfo.pm
#

use strict;
use warnings;
use BeginPerlBioinfo;

my $genbank_library = 'library.gb';

my $gbfh;
unless(open($gbfh,"$genbank_library")) {
	print "Cannot open GenBank library $genbank_library\n";
	exit;
}

my %count = ();
my $size = 4;

while(my $record = get_next_record($gbfh)) {

	# Extract DNA sequence
	my (undef, $dna) = get_annotation_and_dna($record);
	
	# Get count of tetramers
	mercount($size, $dna, \%count);
}

close($gbfh);

#
# Sort the keys by the count, and output results
#
my @sortedkeys = sort {$count{$b} <=> $count{$a}} keys %count;

foreach my $key (@sortedkeys) {
	print "$key ", $count{$key}, "\n";
}

exit;

################################################################################
# Subroutines
################################################################################
# mercount
#  - count all "mers" - subsequences of specified size - in a sequence
#
sub mercount {
	my($size, $seq, $count) = @_;
	# $count is a pointer to a hash

	# Iterate through each subsequence
	for(my $i=0 ; $i <= length($seq)-$size ; ++$i) {

		my $mer = substr($seq, $i, $size);

		if(defined $$count{$mer}) {
			$$count{$mer}++;
		}else{
			$$count{$mer} = 1;
		}
	}
}
